
diffmean_helper <- function(data, treatment, item){

  treat_filter <- data[,treatment] == 1 
  control_filter <- data[,treatment] == 0 
  treat <- subset(data, treat_filter) 
  control <- subset(data, control_filter)
  point <- mean(treat[,item], na.rm=T) - mean(control[,item], na.rm=T)
  n1 <- sum(!is.na(treat[,item]))
  n0 <- sum(!is.na(control[,item]))
  se <- sqrt(var(treat[,item], na.rm=T)/n1 + var(control[,item], na.rm=T)/n0)
  #upper <- point + qnorm(0.975)*se
  #lower <- point + qnorm(0.025)*se
  return(c(point, se))
}

diffmean_simple <- function(data, treatment, item){

  out <- diffmean_helper(data, treatment, item)
  point <- out[1]; se <- out[2]
  upper <- point + qnorm(0.975)*se
  lower <- point + qnorm(0.025)*se
  p <- 2*(1-pnorm(abs(point/se)))
  return(c(lower, point, upper, p, se))
}

diffmean_block <- function(data, treatment, item){
  data_ls <- split(data, data$block_rand)
  group_est <- as.data.frame(t(sapply(data_ls, diffmean_helper, treatment, item)))
  colnames(group_est) <- c("ate", "se")
  block_n <- as.numeric(table(data$block_rand)) # # of obs in each block
  group_est$n <- block_n
  N <- sum(group_est$n)

  # weighted sum for point estimate
  point <- sum(group_est$ate * group_est$n / N)
  # weighted sum for se
  se <- sqrt(sum(group_est$se^2 * (group_est$n/N)^2))
  
  upper <- point + qnorm(0.975)*se
  lower <- point + qnorm(0.025)*se
  p <- 2*(1-pnorm(abs(point/se)))
  return(c(lower, point, upper, p, se))
  
}

